Is efficiency of classical simulations of quantum dynamics related to integrability? 
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Efficiency of time-evoiution of quantum observabies, and thermai states of quenched hamiltonians, 
is studied using time- dependent density matrix renormalization group method in a family of generic 
quantum spin chains which undergo a transition from integrable to non-integrable - quantum chaotic 
case as control parameters are varied. Quantum states (observabies) are represented in terms of 
matrix-product-operators with rank D E (t), such that evolution of a long chain is accurate within 
fidelity error e up to time t. We find that the rank generally increases exponentially D e (t) oc 
exp(constt), unless the system is integrable in which case we find polynomial increase. 
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In the theory of classical dynamical systems there is 
a fundamental difference between integrable and chaotic 
systems. Chaotic systems, having positive algorithmic 
complexity, unlike the integrable ones, cannot be simu- 
lated for arbitrary times with a finite amount of informa- 
tion about their initial states. Computational complexity 
of individual chaotic trajectories is linear in time, how- 
ever, if one wants to describe statistical states (phase 
space distributions) or observabies of chaotic classical 
systems, up to time t, exponential amount of computa- 
tional resources N(t) ~ exp(ht) is needed, where h is the 
Kolmogorov's dynamical entropy related to exponential 
sensitivity to initial conditions. For example, one needs 
to expand the solution of the Liouville equation into the 
lowest N(t) Fourier modes. 

How difficult is it to simulate isolated and bounded 
quantum systems of many interacting particles us- 
ing classical resources? In analogy with the classical 
(chaotic) case, we might expect that the best classical 
simulation of typical quantum systems (in thermody- 
namic limit (TL)) is exponentially hard. Even though 
there is no exponential sensitivity to initial conditions in 
quantum mechanics, there is a tensor-product structure 
of the many-body quantum state space which makes its 
dimension to scale exponentially with the number of par- 
ticles, as opposed to linear scaling in the classical case. 
Furthermore, due to intricate quantum correlations (en- 
tanglement) generic quantum time evolution cannot be 
reduced to (efficient) classical computation in terms of 
non-entangled (classical like) states. However, it is not 
known what amount and form of quantum entanglement 
is needed in order to prevent efficient classical simulation. 

Recently, a family of numerical methods for the simu- 
lation of interacting many-body systems has been devel- 
oped yj which is usually referred to as time-dependent 
density- matrix-renormalization group (t-DMRG), and 
which has been shown to often provide an efficient clas- 
sical simulation of certain interacting quantum systems. 
Simulations of locally interacting one-dimensional quan- 
tum lattices were actually shown rigorously to be efficient 



in the number n of particles Q (i.e., computation time 
and memory resources scale as polynomial functions of n 
at fixed t, or up to t = O(logn)), whereas the scaling of 
computation time and memory with physical time t (in 
TL n = oo), later on referred to as time efficiency, has 
not been systematically studied. t-DMRG was shown to 
be time efficient 4 3j only in rather special cases of exactly 
solvable dynamics (generated with XY spin chain Hamil- 
tonian) and/or for particular choices of initial states, ly- 
ing either in low-energy-sectors or in low dimensional 
invariant subspaces. However, for applications in non- 
equilibrium statistical mechanics and condensed matter 
theory, e.g. in transport phenomena, it is of primary 
importance to understand long-time dynamics of generic 
interacting quantum systems |^|. 

In this Letter we address the question of time effi- 
ciency implementing up-to-date version of t-DMRG for 
a family of Ising spin- 1/2 chains in arbitrary oriented 
magnetic field, which undergoes a transition from in- 
tegrable (transverse Ising) to non-integrable quantum 
chaotic regime as the magnetic field is varied. We fo- 
cus on evolution of density operators of mixed states, 
starting from a thermal state of a quenched hamilto- 
nian, and evolution of local or extensive initial observ- 
abies in Heisenberg picture. Note that time evolution of 
pure states is often ill defined in TL Q. As a quanti- 
tative measure of time efficiency we define and compute 
the minimal dimension D e (t) of matrix product opera- 
tor (MPO) representation of quantum states/observables 
which describes time evolution up to time t within fi- 
delity 1 — 0(e). Our central result states that in generic 
non-integrable cases computation resources grow expo- 
nentially D e (t) oc exp(/i q i), except in the integrable case 
of transverse Ising chain, where the growth is typically 
linear D e (t) oc t. Constant h q , asymptotically indepen- 
dent of n, depends only on the evolution (hamiltonian) 
and not on the details of the initial state/observable or 
error measures, and can be interpreted as a kind of quan- 
tum dynamical entropy. We conjecture that integrability 
(solvability) of Id interacting quantum systems is in one- 
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to-one correspondence with the efficiency of their classi- 
cal simulability. 

We also studied time efficiency of simulation of pure 
states in Schrodinger picture, for which many examples 
of efficient applications exist, however all for initial states 
of rather particular structure typically corresponding to 
low energy sectors. Treating other, typical states, e.g. 
eigenstates of unrelated Hamiltonians, linear combina- 
tions of highly excited states, or states chosen randomly 
in the many-particle Hilbert space, we found that, irre- 
spectively of integrability of dynamics, t-DMRG is not 
time-efficient, i.e. D e (t) grows exponentially even in the 
integrable case of transverse field (consistently with a lin- 
ear growth of entanglement entropy In view of this 
fact, our finding that t-DMRG can be time-efficient for 
integrable systems when implemented for time-evolved 
operators or high-temperature thermal states, provides a 
new paradigm fora successful application of t-DMRG. 

Let us briefly review t-DMRG for evolution of density 
matrices and operators Q which generalizes t-DMRG for 
pure states Q. One defines a superket corresponding to 
an operator O expanded in the computational basis of 
products of local operators. Concretely, for a chain of n 
qubits we use a basis of 4™ Pauli operators a s ° ® • • • ® 
(J s„_i^ s . g {0,x, y, z} and <7° = 1. The key idea of 
t-DMRG is to represent any operator in a matrix prod- 
uct form, O = Y, s j tr (^o° ' ' ' A n-i ) °" S ° ® ' ' ' ® . in 
terms of An matrices Aj 3 of fixed dimension D. The num- 
ber of parameters in the MPO representation is AnD 2 and 
for sufficiently large D it can describe any operator. In 
fact, the minimal D required equals to the maximal rank 
of the reduced super-density-matrix over all bi-partitions 
of the chain. The advantage of MPO representation lies 
in the fact that doing an elementary local one or two 
qubit unitary transformation O' = WOU can be done 
locally, affecting only a pair of neighboring matrices A*? . 

In order to study the role of integrability on the effi- 
ciency of t-DMRG we take antiferromagnetic Ising chain 
in a general homogeneous magnetic field, 

n-2 n— 1 

H(h*, h') = ofcf +1 + 5>V* + h'ofr (1) 

3=0 j=0 

where <r| = 1®-? ® a s <g> We will analyze 

evolution for two different magnetic field values: (i) an 
integrable (regular) case Hr — H(Q, 2) with transverse 
magnetic field and (ii) non-integrable (quantum chaotic) 
case He = H(l, 1) with tilted magnetic field. Particular 
value of h z = 2 in the case of -Hr plays no role. -Hr 
can be solved by Jordan- Wigner transformation which 
maps Hr to a system of noninteracting fermions. To 
confirm that He, and -Hr, indeed represent generic quan- 
tum chaotic, and regular, system, respectively, we calcu- 
lated level spacing distribution (LSD) of their spectra 
(shown in Fig. ^| . LSD is a standard indicator of quan- 
tum chaos 13 . It displays characteristic level repulsion for 



strongly non-integrable quantum systems, whereas for in- 
tegrable systems there is no repulsion due to existence of 
conservation laws and quantum numbers. 




FIG. 1: Nearest neighbor LSD for H c (left) and H R (right) for 
n — 12. Dashed curves are p(s) = S7r/2exp (-Ti-V/4) (left) 
and p(s) = exp( — s) (right), typical for chaotic and regular 
systems, respectively Q. Eigenenergies £ [—9,9] were used 
and statistics for even and odd parity states were combined. 

Evolution by t-DMRG proceeds by splitting hamil- 
tonian (JIJ into even and odd terms, H = H c + H Q , 
such that terms within H e or H a commute between each 
other. An approximate propagator for short time-step 
is then written using Trotter-Suzuki formula as U(St) — 
c -iHc8t/2 e -iH 6t e -iHc8t/2 ^ wnere eac h Q f ^ e three terms 

can be written as a series of commuting one and two 
qubit operations. There are two sources of errors in t- 
DMRG scheme. One is Trotter error scaling as cx (5t) 3 
per time step, or oc (St) 2 in total, and the other, usually 
dominating one, is due to truncation. 

The truncation error arrises because after performing 
two qubit transformation on MPO the required dimen- 
sion of the new matrices increases to AD. In order to pre- 
vent the exponential growth of D with time we truncate 
the resulting matrices back to dimension D Q. Trun- 
cation after application of a single gate Ui introduces a 
norm error n(Ui) equal to the sum of squares of discarded 
singular values. As an estimate for the total truncation 
error jftot (t) at time t we will use a sum of all truncation 
errors rj(Ui) for two qubit gates Ui applied upto time t, 
U(t) — JT, U (the number of such gates scales as ~ t/St). 
If A|(J7j), k = 0, . . . AD — 1, denote decreasingly ordered 
eigenvalues of the reduced super-density-matrix after the 
application of a gate Ui, then 77(Z7j) and 7y to t are given by 

4D-1 

t](Ui) - ]T X 2 (U), Vtot (t) =$>(^)- (2) 

j=D i 

Simple perturbation argument shows that for small time 
step St, single gate truncation error scale as r)(Ui) oc 
(St) 2 , so the total error r] to t(t) oc St. We use the same 
time step St = 0.01 in all our simulations. One may hope 
that Tftot (t) gives a good measure of fidelity 

|tr{0 M po(t)0 exact 
U |tr{0^ PO W}l|tr{O c 2 xact (t)}|' 1 ' 

where Ompo(^) is an operator obtained from the initial 
O with t-DMRG evolution with a given fixed D, while 
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O cxact (t) = W(t)OU{t) is obtained with an exact evo- 
lution. Indeed, by comparing to exact numerical sim- 
ulations of small systems of size n — 6, 8, 10 and sev- 
eral different D we find quite generally (see Fig. [21 for 
an example) that up to good numerical approximation 
1 — F{t) « critot(t) / 5t, where c is some numerical con- 
stant of order 1 which does not depend on 5t, D or n. 



does not depend on e, properties of O(0) or n, for big 




FIG. 2: Fidelity J3J of t-DMRG evolution (full curves) and 
scaled truncation errors crjtot(t)/8t with c = 0.5 (dashed 



curves), for O(0) = cI/ 2 , hamiltonian Hb. and n = 10. Differ- 
ent sets of curves are for D = 10,20,30,40 (top to bottom). 
Chain line marks the threshold where the truncation error 
rjtot{t) = 10~ 4 (indicated by stars for different D's). 

The central quantity we are going to study is D e (t) 
which is the minimal dimension D of matrices A?* in 
order for the total truncation error ijtot {t) to be less than 
some error tolerance e, or fidelity 10 to be bigger than 1 — 
(c/St)e, for evolution to time t. We use e = 10~ 4 for local 
and extensive operators and e = 10 -6 for thermal states. 
The central question is: does D e (t) grow exponentially 
or polynomially with t? If it grows polynomially we can 
say that t-DMRG is time efficient. 

Let us first study the case where the initial operator is a 
local operator in the center of the lattice O(0) = cr^ 2 , s <E 
{x, y,z}. In the integrable case time evolution 0(t) can 
be computed exactly in terms of Jordan- Wigner trans- 
formation and Toeplitz determinants || , however for ini- 
tial operators with infinite index [1(1 ]. like e.g for o~*'J 2 , 
n — > oo, the evolution is rather complex and the effective 
number of terms (Pauli group elements) needed to span 
Oft) grows exponentially in t. In spite of that, our nu- 
merical simulations shown in Fig. [31 strongly suggest the 
linear growth D e (t) ~ t for initial operators with infinite 
index. Quite interestingly, for initial operators with finite 
index, D e (t) saturates to a finite value, for example llj 



D e (oo) = 4 for ct^ /2 , or D e (oo) = 16 for a^ l/2 _ 1 a] 



Va- 



in 



non-integrable cases the rank has been found to grow ex- 
ponentially, D £ (t) ~ exp(h q t) with exponent h q which 



For H = H c we find h a = 1.10 
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FIG. 3: D e (i) for local initial operators. We consider three 
cases O(0) = a^J^ (empty circles, squares and triangles), for 
non-integrable evolution He, and four cases, O(0) = &*'J 2 (full 
squares, diamonds), cr n/2-i <7 n/2 (^ uu triangles) with infinite 
index, and O(0) = o _ ^/2-i cr n/2 (f uu circles) with index 2, 
for integrable evolution Hr. Full line in the inset illustrates 
exponential growth oc e l u in the non-integrable case. Full 
squares and diamonds are for n = 40, otherwise n = 20. 



In physics it is often useful to consider extensive ob- 
servables, for instance translational sums of local oper- 
ators, e.g. the hamiltonian H or the total magnetiza- 
tion M s — Y^j=Q a j- A s opposed to local operators, 
extensive initial operators, interpreted as W-like super- 
states, contain some long-range entanglement so one may 
expect that t-DMRG should be somewhat less efficient 
than for local operators. Indeed, in the integrable case 
we find for extensive operators with finite index that 
D e (t) does no longer saturate but now grows linearly, 
D e (t) ~ t, whereas for extensive operators with infinite 
index the growth may be even somewhat faster, most 
likely quadratic D e (t) ~ t 2 but clearly slower than ex- 
ponential. In the non-integrable case, we again find ex- 
ponential growth D e (t) ~ exp(/i q t) with the same expo- 
nent h q as for local initial observables. The results are 
summarized in Fig. Note that for local as well as for 
extensive observables 7? to t (t) asymptotically does not de- 
pend on n. Therefore the results shown in Figs. I.'il II for 
which convergence with n has been reached, are already 
representative of TL. 

In the last set of numerical experiments we consider 
time efficiency of the evolution of a thermal initial state 
0(t) = Z^ 1 exp(— (3Hq) under a sudden change of the 
hamiltonian at t = 0, namely H(t < 0) = Hq = 
H(0, l),H(t > 0) = Hi. Again, we treat two situa- 
tions: in the first case we consider change after which 
the hamiltonian remains integrable, Hi = H(0, 2) = Hr , 
while in the other case the change breaks integrability 
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FIG. 4: D ( (t) for extensive initial operators. For both hamil- 
tonians He, Hb we take 0(0) = ^ .a* (empty, full squares) 
with infinite index, and 0(0) = H(0, 1) (empty, full cir- 
cles) with index 1. For Hb we also show case 0(0) = 
y\ cr|<j^ +1 + o"Jo"J + i (full diamonds) with index 1 and 2. In 
the semi- log inset we illustrate exponential increase oc e * 
(full straight line) for He and polynomial ~ t 2 (full curve) for 
Hb- For full circles n = 64, otherwise n = 32. 



longer entanglement range as we increase the power p. 
These results are summarized in Fig. [5] In contrast to 
local and W-like observables, the total truncation error 
Vtot{t) is f° r thermal states proportional to n. Therefore, 
the fidelity at fixed t and D of t-DMRG simulation of 
thermal states decreases in TL. 

In conclusion, we have presented numerical experi- 
ments suggesting that the scaling of classical computa- 
tion resources in t-DMRG simulations of quantum Id lat- 
tices with local interaction may sensitively depend on the 
integrability of the hamiltonian, and on whether we prop- 
agate pure states or mixed states/observables. For the 
latter we find universal exponential growth of the mini- 
mal rank of the matrix product representation in physical 
time, unless wc propagate by an integrable hamiltonian 
from the initial state/observable which can be related to 
(sums of) local operators, in which case the growth is 
polynomial, or even saturates for a specific class of ini- 
tial operators. We acknowledge stimulating discussions 
with J. Eisert, A.J. Daley, and P. Zoller, and support 
by Slovenian Research Agency, programme P 1-0044, and 
grant Jl-7437. 
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FIG. 5: D e (t) for thermal states of H with ft = 0.01 
(/3 = 0.05 in inset), for evolution with He (open symbols) 
and Hb, (full symbols) at n = 40. Solid curves again indicate 
exponential increase oc e . 



of the hamiltonian, Hi = H(l, 1) = He- Initial state is 
prepared from identity super-state using imaginary time 
t-DMRG with the same MPO rank D as it is later used 
for real time dynamics. We find, consistently with previ- 
ous results, that at high temperature (/3 < 1) the rank 
D £ (t) grows very slowly, perhaps slower than linear, in 
the integrable case, and exponentially D e (t) ~ exp(/i q t), 
in the non-integrable case. Interestingly, at lower tem- 
peratures we find exponential growth in both cases, even 
in the integrable one. This is not unreasonable as the 
initial (thermal) state can be expanded in a power series 
in (3 and the higher orders become less local with 
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